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Abstract 

We  present  a  generalized  polynomial  chaos  algorithm  for  the  solution  of  stochastic  elliptic  partial  differential 
equations  suject  to  uncertain  inputs.  In  particular,  we  focus  on  the  solution  of  the  Poisson  equation  with  ran¬ 
dom  diffusivity,  forcing  and  boundary  conditions.  The  stochastic  input  and  solution  are  represented  spectrally 
by  employing  the  orthogonal  polynomial  functionals  from  the  Askey  scheme,  as  a  generalization  of  the  original 
polynomial  chaos  idea  of  Wiener  (1938).  A  Galerkin  projection  in  random  space  is  applied  to  derive  the  equations 
in  the  weak  form.  The  resulting  set  of  deterministic  equations  for  each  random  mode  is  solved  iteratively  by  a 
block  Gauss-Seidel  iteration  technique.  Both  discrete  and  continuous  random  distributions  are  considered,  and 
convergence  is  verified  in  model  problems  and  against  Monte  Garlo  simulations. 
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1  Introduction 


It  has  been  common  practice  in  engineering  to  analyze  systems  based  on  deterministic  mathematical  models  with 
precisely  defined  input  data.  However,  since  such  ideal  situations  are  rarely  encountered  in  practice,  the  need 
to  address  uncertainties  is  now  clearly  recognized,  and  there  has  been  a  growing  interest  in  the  applications  of 
probabilistic  methods  [1,  2,  3,  4]. 

The  probabilistic  methods  in  engineering  can  be  broadly  classified  into  two  major  categories:  methods  using  a 
statistical  approach  and  methods  using  a  non-statistical  approach.  The  statistical  approach  includes  Monte  Carlo 
simulation,  stratified  sampling,  Latin  hypercube  sampling,  etc.  These  methods  involve  sampling  and  estimation  and 
in  most  cases  are  straightforward  to  apply.  However,  since  the  accuracy  of  the  sampling  techniques  depends  on  the 
sample  size,  in  accordance  with  the  ‘weak  law  of  large  number’,  simulations  can  become  prohibitively  expensive, 
especially  for  the  systems  that  are  already  complicated  in  the  deterministic  case.  Thus,  these  methods  are  often  used 
as  the  last  resort  and  research  effort  has  been  made  in  developing  the  non-statistical  methods. 

The  most  popular  non-statistical  method  is  the  perturbation  method,  where  the  random  field  is  expanded  via 
Taylor  series  around  its  mean  and  truncated  at  certain  order.  Typically,  at  most  second-order  expansion  is  employed 
because  the  system  of  equations  becomes  extremely  complicated  beyond  second-order.  This  approach,  also  called 
the  ‘second  moment  analysis’  [5,  6,  7],  has  been  used  extensively  in  various  fields  [8].  An  inherent  limitation  of  the 
perturbation  method  is  that  the  uncertainties  cannot  be  too  large,  i.e.  variances  of  the  random  field  cannot  be  too 
large  compared  with  their  mean  values,  e.g.  typically  less  than  10%  [5].  Also,  higher-order  statistics  are  not  readily 
available  for  the  second  moment  method.  Another  approach  is  the  Neumann  expansion,  which  is  based  on  the  inverse 
of  the  stochastic  operator  in  a  Neumann  series.  This  method  too  is  restricted  to  small  uncertainties  and  attempts 
have  been  made  to  couple  it  with  the  Monte  Carlo  simulation  to  obtain  more  efficient  algorithms  [9,  10]. 

Another  methodology  of  the  non-statistical  type  is  to  ‘discretize’  directly  the  random  field.  Ghanem  &  Spanos 
pioneered  a  polynomial  chaos  expansion  method  and  have  successfully  applied  it  to  various  problems  in  mechanics 
[11].  The  polynomial  chaos  expansion  is  based  on  the  homogeneous  chaos  theory  of  Wiener  [12]  and  is  essentially 
a  spectral  expansion  of  the  random  variables.  It  allows  high-order  representation  and  promises  fast  convergence; 
coupled  with  Karhunen-Loeve  decomposition  for  the  input  and  Galerkin  projection  in  random  space,  it  results  in 
computationally  tractable  algorithms  for  large  engineering  systems  [13].  More  efficient  Monte  Carlo  algorithms  can 
also  be  designed  when  combined  with  the  chaos  expansion  technique  [13,  14].  More  recently,  a  theoretical  framework 
of  discretizing  the  random  field  via  the  finite  element  approach,  i.e.  piecewise  polynomials,  was  proposed  in  [15]. 
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The  classical  polynomial  chaos  expansion  is  based  on  the  Hermite  polynomials  in  terms  of  Gaussian  random 
variables.  Although  in  theory,  it  converges  to  any  functionals  in  the  random  space  [16],  it  achieves  optimal 
convergence  rate  only  for  Gausian  and  near  Gaussian  random  fields  [17],  and  does  not  readily  apply  to  the  random 
fields  with  discrete  distribution.  A  more  general  framework,  called  the  ‘generalized  polynomial  chaos’  or  the  ‘Askey- 
chaos’,  was  proposed  in  [17].  Here  the  polynomials  are  chosen  from  the  hypergeometric  polynomials  of  the  Askey 
scheme  [18],  and  the  underlying  random  variables  are  not  restricted  to  Gaussian  random  variables.  Instead,  the 
type  of  random  variables  are  chosen  according  to  the  stochastic  input  and  the  weighting  function  of  these  random 
variables  determines  the  type  of  orthogonal  polynomials  to  be  used  as  the  basis  in  the  random  space.  The  convergence 
properties  of  different  bases  were  studied  in  [17]  and  exponential  convergence  rate  was  demonstrated  for  model 
problems.  In  [19]  the  generalized  polynomial  chaos  was  applied  to  modeling  uncertainties  in  incompressible  flow  and 
in  [20]  to  flow-structure  interaction  problems. 

The  main  objective  of  this  paper  is  to  give  a  broad  algorithmic  framework  to  solve  stochastic  elliptic  partial 
differential  equations  based  on  the  generalized  polynomial  chaos  expansion.  The  class  of  problems  we  solve  has  the 
form 

J  W  ■  [k{x]  oj)W u{x]  oj)]  =  f  {x]  oj) ,  (x;  w)  G  H  X  n  ,  , 

^  u{x;lo)  =  g{x;Lo),  {x;u!)  €  dD  x  ' 

where  U  is  a  bounded  domain  in  {d  =  1, 2, 3)  and  is  a  probability  space.  /,  g  and  k  are  M-values  functions  on 

D  X  n.  This  can  be  considered  as  a  model  of  steady  state  diffusion  problems  subject  to  internal  (diffusivity  k)  and/or 

external  (source  term  /  and/or  Dirichlet  boundary  condition  g)  uncertainties.  Babuska  was  among  the  first  to  study 

rigorously  existence  of  solutions  of  the  random  Dirichlet  problem  [21].  Becus  &  Gozzarelli  studied  the  existence  and 

properties  of  the  general  solution  to  (1),  see  [22,  23,  24].  Also,  in  [15]  the  problem  subject  to  random  diffusivity 

and/or  random  source  terms  was  studied  and  existence  and  uniqueness  of  the  solution  in  the  finite  element  concept, 

both  in  physical  space  and  random  space,  were  addressed.  From  a  different  perspective,  problems  similar  to  (1)  were 

studied  via  homogenization  methods  to  evaluate  the  ‘effective  diffusivity’  of  random  media  [25] . 

In  this  paper,  we  solve  the  steady  state  diffusion  problem  (1)  by  generalized  polynomial  chaos  expansion,  where 
the  uncertainties  can  be  introduced  through  k,  f,org,  or  some  combinations.  It  is  worth  noting  that  when  both 
K  and  u  are  random,  it  is  not  obvious  how  to  give  a  mathematical  meaning  or  justification  to  the  product  of  two 
stochastic  processes  if  they  are  not  smooth.  However,  the  product  is  well  defined  in  terms  of  the  chaos  expansion  by 
using  the  concept  of  Wick  product  and  Kondratiev  space  [26,  27]. 

In  the  next  section,  we  present  the  concept  of  the  generalized  polynomial  chaos,  and  in  section  3  we  apply  the 
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expansion  to  the  solution  of  the  steady  state  diffusion  problem.  Numerical  results  are  presented  in  section  4,  and  we 
conclude  the  paper  with  a  discussion. 


2  The  Generalized  Polynomial  Chaos 


In  this  section  we  introduce  the  generalized  polynomial  chaos  expansion  along  with  the  Karhunen-Loeve  (KL)  de¬ 
composition,  another  classical  technique  for  representing  random  processes.  The  KL  decomposition  can  be  used  in 
some  cases  to  represent  efficiently  the  known  stochastic  fields,  i.e.,  the  stochastic  inputs. 

2.1  The  Askey  Scheme 


The  generalized  hypergeometric  series  rFs  is  defined  by 


OO 

rFs{ai,---  ,ar;bi,---  ,bs]z)  =  ^ 


(ai)fc  ■  ■  ■  (ar)fc  z’^ 
{bi)k---{bs)k  fc!’ 


(2) 


where  bi  yf  0,  —  1,  —2, . . .  for  t  =  {1, . . . ,  s}  to  ensure  the  denominator  factors  in  the  terms  of  the  series  are  not  zero. 
The  Pochhammer  symbol  (a)„  defined  as 


r  1,  if  n  =  0, 

a(a -I- 1)  •  •  •  (a -I- n  —  1),  if  n  =  1, 2,  3, . . . . 


(3) 


If  one  of  the  numerator  parameters  Uj,  i  =  I, . . . ,  r  is  a  negative  integer,  say  ai  =  — n,  the  hypergeometric  series  (2) 
terminates  at  the  n*^-term  and  becomes  a  hypergeometric  polynomial  in  z. 

The  Askey  scheme,  which  is  represented  as  a  tree  structure  in  figure  I  (following  [28]),  classifies  the  hypergeometric 
orthogonal  polynomials  and  indicates  the  limit  relations  between  them.  The  ‘tree’  starts  with  the  Wilson  polynomials 
and  the  Racah  polynomials  on  the  top.  The  Wilson  polynomials  are  continuous  while  the  Racah  polynomials  are 
discrete.  The  lines  connecting  different  polynomials  denote  the  limit  transition  relationships  between  them;  this 
implies  that  the  polynomials  at  the  lower  end  of  the  lines  can  be  obtained  by  taking  the  limit  of  one  of  the  parameters 
from  their  counterparts  on  the  upper  end.  For  example,  the  limit  relation  between  Jacobi  polynomials  Pn°’^\x)  and 
Hermite  polynomials  iL„(x)  is 


lim 

a— ^OO 


Hnjx) 
2”n!  ’ 


and  between  Meixner  polynomials  M„(x;/3,  c)  and  Charlier  polynomials  C„(a;;a)  is 

lim  M„  (  x-,P,  )  =  Cnix;a). 

/3^oo  \  a  +  p  J 

For  a  detailed  account  of  definitions  and  properties  of  hypergeometric  polynomials,  see  [18];  for  the  limit  relations  of 
Askey  scheme,  see  [29]  and  [28]. 
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Figure  1:  The  Askey  scheme  of  orthogonal  polynomials 
The  orthogonal  polynomials  associated  with  the  generalized  polynomial  chaos,  include: 

•  Hermite,  Laguerre  and  Jacobi  polynomials  for  continuous  distributions; 

•  Charlier,  Meixner,  Krawtchouk  and  Hahn  polynomials  for  discrete  distributions. 

2.2  The  Generalized  Polynomial  Chaos:  Askey-Chaos 

The  original  polynomial  chaos  [12,  30]  employs  the  Hermite  polynomials  in  the  random  space  as  the  trial  basis  to 
expand  the  stochastic  processes.  Cameron  &  Martin  proved  that  such  expansion  converges  to  any  second-order 
processes  in  the  L2  sense  [16].  It  can  be  seen  from  figure  1  that  the  Hermite  polynomials  are  a  subset  of  the  Askey 
scheme.  The  generalized  polynomial  chaos,  or  the  Askey-Chaos,  was  proposed  in  [17,  19,  20]  and  employs  more  types 
of  orthogonal  polynomials  from  the  Askey  scheme.  Convergence  to  second-order  stochastic  processes  can  be  possibly 
obtained  as  a  generalization  of  Cameron-Martin  theorem  [16]. 

The  chaos  expansion  is  essentially  a  representation  of  a  function  /  G  T2(f^)  where  is  the  properly  defined 
probability  space.  We  denote  by  orthogonal  polynomials  from  the  Askey  scheme,  which  form  an 

orthogonal  basis  in  T2(R")- 
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A  general  second-order  random  process  X(lu)  can  be  represented  in  the  form 


X{u})  =  co'I'o 

oo 

OO  ii 

il  —  l  22  =  1 
oo  il  12 

+  EEE  ^212223  ^3 

*1  =  1  *2  =  1  *3  =  1 

+  •••,  (4) 

where  denotes  the  generalized  polynomial  chaos  of  order  n  in  terms  of  the  multi-dimensional  random 

variables  ^  =  (^i, . . . , . . . ).  Note  that  this  is  an  infinite  summation  in  the  infinite  dimensional  space  of  The  ex¬ 
pansion  bases  {ihr*}  are  multi-dimensional  hypergeometric  polynomials  defined  as  tensor-products  of  the  correspond¬ 
ing  one-dimensional  polynomials  bases  X  be  the  space  of  index  sequences  (ai,  02,  •  •  ■ ,  ocn,  •  ■  • )  €  Nq 

and  n  :=  otk-  Then 

n 

4'**(a,6,---,er*)=  (5) 

fc=l 

For  notational  and  computational  convenience,  equation  (4)  is  often  rewritten,  according  to  some  numbering  scheme, 
in  the  form  with  only  one  index  as 

OO 

j=0 

where  there  is  a  one-to-one  correspondence  between  the  coefficients  and  basis  functions  in  (4)  and  (6). 

The  family  {4)„}  is  an  orthogonal  basis  in  L2{^)  with  orthogonality  relation 

<  >=<  >  5ij,  (7) 

where  5ij  is  the  Kronecker  delta  and  <  •,  •  >  denotes  the  ensemble  average  which  is  the  inner  product  in  the  Hilbert 
space  of  the  variables  ^ 


<  >=  I  fiOgimmt  (8) 

or 

</(l)5(0>=E/(^)5(0Vb(4)  (9) 

€ 

in  the  discrete  case.  Here  kF(^)  is  the  weighting  function  corresponding  to  the  Askey  polynomials  chaos  basis  {4)^}. 


6 


In  the  original  polynomial  chaos,  {‘hn}  are  the  Hermite  polynomials  and  ^  are  the  Gaussian  random  variables. 
In  the  Askey-chaos  expansion,  the  orthogonal  polynomials  {d>„}  are  not  restricted  to  Hermite  polynomials,  instead 
they  are  determined  by  the  weighting  function  of  the  corresponding  random  variables  which  are  not  necessarily 
Gaussian  variables.  The  corresponding  type  of  polynomials  and  their  associated  random  variables  are  listed  in  table 
1.  Since  each  type  of  polynomials  from  the  Askey  scheme  form  a  complete  basis  in  the  Hilbert  space  determined  by 


Random  variables  ^ 

Orthogonal  polynomials  {$„} 

Support 

Continuous 

Gaussian 

Hermite 

(—00, 00) 

Gamma 

Laguerre 

[0,oo) 

Beta 

Jacobi 

[a,b] 

Uniform 

Legendre 

[a,b] 

Discrete 

Poisson 

Charlier 

{0,1,2,...} 

Binomial 

Krawtchouk 

{0,1,. ..,1V} 

Negative  Binomial 

Meixner 

{0,1,2,...} 

Hypergeometric 

Hahn 

{0,1,. ..,1V} 

Table  1:  Correspondence  of  the  orthogonal  polynomials  and  random  variables  for  different  Askey-chaos  {N  >  0  is  a 
finite  integer). 


their  corresponding  support,  we  can  expect  each  type  of  Askey-chaos  to  converge  to  any  L2  functional  in  the  L2  sense 
in  the  corresponding  Hilbert  functional  space  as  a  generalized  result  of  the  Cameron-Martin  theorem  [16,  31].  Each 
type  of  orthogonal  polynomials  from  the  Askey-chaos  has  weighting  functions  of  the  same  form  as  the  probability 
function  of  its  associated  random  variables  as  shown  in  table  1. 

The  original  polynomial  chaos,  the  Hermite-chaos,  is  a  subset  of  the  continuous  chaos  and  has  been  studied  in 
the  literature  extensively.  Here  we  show  the  Laguerre-chaos  as  another  example  of  the  continuous  chaos.  In  one 
dimension,  the  nth-order  Laguerre  polynomials  are  defined  as 

(C)  =  (I) ”  (e-«r+“)  ,  e  >  0,  (10) 


where  parameter  a  >  —1  is  a  real  positive  number.  The  weight  function  in  the  orthogonality  relation  (7)  is 


W{i-a)  = 


T{a  +  1)' 


(11) 


It  can  be  seen  that  this  is  the  same  as  the  probability  density  function  (PDF)  of  a  Gamma  random  variable.  The 
first  few  members  of  the  one-dimensional  Laguerre-chaos  are: 


$0  =  1,  $i=C-(a+l), 


$2  —  -  ^{a  -I-  2)  -I-  -{a  +  l)(a  -I-  2), 


(12) 


As  an  example  of  the  discrete  chaos,  the  Charlier-chaos  employs  the  Charlier  orthogonal  polynomials  as  the  basis 
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in  the  random  space.  In  one  dimension,  the  nth-order  Charlier  polynomial  can  be  written  as 

=  e  =  0,l,2,...  (13) 

r^O  ^ 

where  A  >  0  and  —  1)  •  •  •  (^  —  r  -I-  1).  The  weighting  function  of  Charlier  polynomials  is  the  probability 

function  of  the  Poisson  distribution  with  mean  A 

m5;A)  =  e-^^.  (14) 

The  first  few  members  of  the  one-dimensional  Charlier-chaos  are: 

$0  =  1,  ^i=e-l,  $2  =e(?-l)-2Ae  +  A",  ...  (15) 

2.3  The  Karhunen-Loeve  Decomposition 

The  Karhunen-Loeve  (KL)  expansion  is  another  way  of  representing  a  random  process  [32].  It  is  based  on  the 
spectral  expansion  of  the  correlation  function  of  the  process.  It  is  particularly  useful  for  the  generalized  polynomial 
chaos  expansion  as  it  provides  a  means  of  reducing  the  dimensionality  of  the  random  space.  Let  us  denote  the 
process  by  /i(x;  w)  and  its  correlation  function  by  Rhh{^,y),  where  x  and  y  are  the  spatial  or  temporal  coordinates. 
By  definition,  the  correlation  function  is  real,  symmetric,  and  positive  definite.  All  eigenfunctions  are  mutually 
orthogonal  and  form  a  complete  set  spanning  the  function  space  to  which  /i(x;  u)  belongs.  The  KL  expansion  then 
takes  the  following  form: 

OO 

/i(x;  oj)  =  /i(x)  -I-  E  \/^</>i(x)C»(w),  (16) 

i=l 

where  /i(x)  denotes  the  mean  of  the  random  process,  and  ^i(w)  forms  a  set  of  uncorrelated  random  variables.  Also, 
4>i{x)  and  Xi  are  the  eigenfunctions  and  eigenvalues  of  the  correlation  function,  respectively,  i.e., 

J  Rhh{^,y)(pi{y)dy  =  Xi(l),{x).  (17) 

Among  other  possible  decompositions  of  a  random  process,  the  KL  expansion  is  optimal  in  the  sense  that  the 
mean-square  error  of  the  finite  representation  of  the  process  is  minimized.  Its  use,  however,  is  limited  as  the  correlation 
function  of  the  solution  process  is  often  not  known  a  priori.  Nevertheless,  the  KL  expansion  provides  an  effective 
means  of  representing  the  input  random  processes  when  the  correlation  structure  is  known. 


3  Galerkin  Projection 


In  this  section  we  present  the  detailed  algorithm  for  the  application  of  the  generalized  polynomial  chaos  expansion 
to  equation  (1).  By  applying  the  chaos  expansion,  we  expand  the  variables  as 

M  M  M 

k{x-,uj)  =  u{x-,uj)  =  '^u^{x)<P^{0,  /(a;;w)  =  (18) 

i— 0  i—0  i— 0 

where  we  have  replaced  the  infinite  summation  of  ^  in  infinite  dimensions  in  equation  (6)  by  a  truncated  finite-term 
summation  of  {<&}  in  the  finite  dimensions  of  ^  The  dimensionality  n  of  ^  is  determined  by  the 

random  inputs.  The  random  parameter  u  is  absorbed  into  the  polynomial  basis  d>(^),  thus  the  expansion  coefficients 
ki  and  ut  are  deterministic.  By  substituting  the  expansion  into  governing  equation  (1),  we  obtain 

M  /  M  \  M 

V-  (19) 

_i=o  \j=o  J  J  i=0 

Upon  simplification,  it  can  be  written  as 

MM  M 

EE[  Ki{x)\/'^Uj{x)  +  \7Ki{x)  ■  \/Uj{x)]  ^  fi{x)^i.  (20) 

i=0  j=0  i=0 

A  Galerkin  projection  of  the  above  equation  onto  each  polynomial  basis  {d>i}  is  then  conducted  in  order  to  ensure 

that  the  error  is  orthogonal  to  the  functional  space  spanned  by  the  finite-dimensional  basis  {d>i}.  By  projecting  with 
<I)fc  for  each  k  =  {0, . . . ,  M}  and  employing  the  orthogonality  relation  (7),  we  obtain  for  each  fc  =  0, . . . ,  M, 

M  M 

E  E  [«i(a;)V^Uj(a;)  -k  VKi{x)  ■  Vuj{x)]  Ciju  =  fk{x)  <  >,  (21) 

i—0  j—0 

where  e^-fc  =<  ‘hi'hj'hfc  >•  Together  with  <  <I)|  >,  the  coefficients  e^-fc  can  be  evaluated  analytically  from  the 
definition  of  <i>i.  By  defining 

M  M 

^3k{x)  =  '^Ki{x)eijk,  hjk{x)  =  Ki{x)eijk  =  Vbjk{x), 

i—0  i—0 

we  can  rewrite  the  above  equation  as 

M 

[bjk{x)\7'^Uj{x)  +  hjk{x)  ■  \/uj{x)]  =  fk{x)  <  >,  Vfc  G  [0,  M].  (22) 

j=o 

Equation  (22)  is  a  set  of  (M  -|-  1)  coupled  elliptic  partial  differential  equations.  These  equations  are  deterministic 
and  can  be  solved  by  any  conventional  method,  e.g.  finite  elements.  In  this  paper  we  employ  the  spectral//ip  element 
method  [33].  The  total  number  of  equations  (M  -|-  1)  is  determined  by  the  dimensionality  of  the  chaos  expansion  (n) 
and  the  highest  order  (p)  of  the  polynomials  {d>},  where 

(M -I- 1)  =  (n -l-p)!/(n!p!).  (23) 
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While  it  is  possible  to  solve  equation  (22)  via  a  direct  solver,  we  choose  to  use  an  iterative  method  to  take  advantage 
of  the  diagonal  dominance  of  the  block  matrix  B  =  {bjk\-  In  particular,  we  employ  the  block  Gauss-Seidel  iteration 
in  the  following  form:  for  all  A:  =  0, . . . ,  M, 


bkk{x)V^ul~^  (x)  +  hkkix)  ■  (x)  =  fk(x)  <  > 

k-l 

-  +  hjk{x)  ■  Vu”+i(x)] 

M 

~  X!  [bjk{x)\7‘^u]{x)  +  hjk{x)  ■  Vu"(x)]  , 

j=k+l 

where  the  superscript  n  denotes  the  iteration  number.  The  convergence  criterion  is  defined  as 

\K+\x)-uUx) 


\u\{x)  —  vPj.{x) 


<  e, 


Vfc  e  [0,M], 


(24) 


(25) 


^kK-^)  '^k\ 

where  e  is  a  small  positive  number  and  different  types  of  norm  ||  •  ||  can  be  used.  In  this  paper  the  Loo  norm  is  used  and 
e  is  set  to  be  10“®  ~  10“^.  For  all  the  results  we  present  here,  the  block  Gauss-Seidel  iteration  normally  converges 
within  about  10  steps.  A  similar  iteration  technique  was  used  in  [34]  for  stochastic  modeling  of  elasto-plastic  body 
problems  with  the  Hermite-chaos  and  fast  convergence  was  reported  too. 


4  Numerical  Results 


In  this  section  we  present  numerical  results  of  the  proposed  generalized  polynomial  chaos  expansion  to  stochastic 
diffusion  problem.  We  first  consider  an  one-dimensional  model  problem  where  the  exact  solution  is  available;  then 
a  more  complicated  two-dimensional  problem  where  we  use  Monte  Garlo  simulation  to  validate  the  chaos  solution. 
Among  the  types  of  chaos  expansions  listed  in  table  1,  we  choose  two  continuous  chaos:  Hermite-chaos  and  Jacobi- 
chaos;  and  two  discrete  chaos:  Gharlier-chaos  and  Krawtchouk-chaos  for  demonstration  purposes.  Finally,  we  solve 
the  random  heat  conduction  problem  in  a  grooved  channel  as  an  example  of  a  more  practical  application. 

4.1  One-Dimensional  Model  Problem 


Gonsider  the  following  problem 


with  boundary  conditions 


d 

dx 


k{x;  to)  —  {x;  to) 
dx 


=  0,  xe[0,l]. 


u(0;uj)  =  0,  m(1;w)  =  1. 


The  random  diffusivity  has  the  form 


k(x;  uj)  =  1  +  e(uj)x, 


(26) 


(27) 
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where  e(w)  is  a  random  variable,  and  k{x',uj)  >  0.  The  exact  solution  to  this  problem  is 


,  In  [1  +  e{oj)x]  /  In  [1  +  e(w)] ,  for  e(w)  ^  0; 

Ue(x,uj)  -  foreH  =  0. 

The  ‘mean-square’  error  of  the  numerical  solution  from  the  generalized  chaos  expansion  Up{x,uj)  is  computed 


(28) 


62 


(x)  =  [Up{x,  Lo) 


where  E  denotes  the  ‘expectation’  operator  and  p  is  the  order  of  the  chaos  expansion.  Specifically,  we  examine 
the  ‘mean-square’  convergence  (L2  convergence  in  random  space)  of  the  Loo  norm  (in  physical  space)  of  62  (x)  as  p 
increases. 

4.1.1  Jacobi-chaos  and  Beta  Distribution 


We  assume  e{uj)  in  equation  (27)  is  a  beta  random  variable,  i.e.  its  probability  density  function  (PDF)  has  the  form 

(l-e)“(l-f  e)^ 


fie;a,P)  = 


2“+/5+iB(a-hl,/3-f  1)^ 


G  [-1,1],  a,p>  -1, 


(29) 


where  B{a,  (3)  is  the  Beta  function  defined  as  B{p,  q)  =  T(p)T{q) /T{p+q).  The  corresponding  generalized  polynomial 
chaos,  according  to  table  1,  is  the  Jacobi-chaos.  An  important  special  case  is  when  a  =  (3  =  0,  then  e(w)  becomes 
an  uniform  random  variable  and  the  corresponding  chaos  becomes  the  Legendre-chaos  (see  table  1). 

In  figure  2  the  mean-square  convergence  of  the  Jacobi-chaos  solution  is  shown  with  different  standard  deviation 
cr  of  the  input  as  a  measure  of  the  magnitude  of  the  input  uncertainty.  It  can  be  seen  on  the  semi-log  scale  that  the 
Jacobi-chaos  solution,  including  the  Legendre-chaos  for  uniform  random  variables,  converges  exponentially  fast  as  the 
expansion  order  p  increases.  The  exponential  convergence  rate  is  retained  for  large  input  uncertainty  such  as  cr  =  0.9, 
which  is  close  to  the  limit  of  the  existence  of  the  solution  (cr  <  1).  This  is  in  contrast  to  the  perturbation-based 
method  which  normally  works  for  a  <  0.1. 


4.1.2  Hermite-chaos  and  Gaussian  Distribution 


We  now  assume  e(u;)  in  equation  (27)  is  a  Gaussian  random  variable  with  probability  density  function 

^  eG  (-00,00).  (30) 

The  corresponding  generalized  polynomial  chaos  is  the  Hermite-chaos  (table  1). 

While  the  random  input  has  infinite  support  and  rigorous  analysis  of  the  existence  and  uniqueness  of  the  solution 
is  lacking  to  ensure  k{x,uj)  >  0  in  equation  (27)  in  some  stochastic  sense,  it  is  intuitive  to  assume  that  the  solution 
exists  for  random  input  with  small  deviation  cr.  In  this  paper,  we  assume  cr  =  0.1  and  the  mean-square  convergence 
of  the  Hermite-chaos  solution  is  shown  in  figure  3.  Again,  exponential  convergence  rate  is  achieved. 
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Figure  2:  Convergence  of  Jacobi-chaos  for  the  one-dimensional  model  problem. 


Figure  3:  Convergence  of  Hermite-chaos  for  the  one-dimensional  model  problem. 
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4.1.3  Charlier-chaos  and  Poisson  Distribution 


We  now  assume  e(u;)  in  equation  (27)  is  a  discrete  random  variable  with  Poisson  distribution 

/(e;A)  =  e-^^,  e  =  0, 1, 2, . . . ,  A  >  0.  (31) 

The  corresponding  generalized  polynomial  chaos  is  the  Charlier-chaos  (table  1).  Again  we  assume  relatively  small 
deviation  ct  =  0.1  to  ensure  the  existence  of  the  solution.  The  exponential  convergence  of  the  Charlier-chaos  expansion 
is  shown  in  figure  4  for  two  different  values  of  the  parameter  A. 


2 
a> 

_r 

1  2  3  4  5 

P 

Figure  4:  Convergence  of  Charlier-chaos  for  the  one-dimensional  model  problem. 

4.1.4  Krawtchouk-chaos  and  Binomial  Distribution 

In  this  section  e{uj)  in  equation  (27)  is  assumed  to  be  a  discrete  random  variable  with  binomial  distribution 

/(e;p,lV)  =  ^^^g"(l  -  g)^"%  0  <  g  <  1,  e  =  0, 1, . . . ,  W  (32) 

The  corresponding  generalized  polynomial  chaos  is  the  Krawtchouk-chaos  (table  1).  Exponential  convergence  of  the 
Krawtchouk-chaos  expansion  can  be  seen  in  figure  5  with  different  values  of  the  parameters  (N,  q) . 

4.2  Two-Dimensional  Model  Problem 

In  this  section  we  consider  the  two-dimensional  problem 

V  ■  [K{x,y;uj)Vu{x,y;uj)]  =  f{x,y;uj),  (x,  y)  G  [-1, 1]  x  [-1, 1]  (33) 

with  boundary  conditions 


u{-l,y,uj)  =  1, 


1^(1, y;w)  =  0, 
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u{x,  — 1;  w)  =  0, 


^(x,  l;w)  =  0. 
dy 


Figure  5:  Convergence  of  Krawtchouk-chaos  for  the  one-dimensional  model  problem. 

The  diffusivity  K{x,y;u))  and  source  term  f{x,y;uj)  are  stochastic  processes  with  certain  distribution  and  given 
correlation  function  C{xi,yi;x2,y2)-  The  mean  fields  are:  K{x,y;uj)  =  1  and  f{x,y;cij)  =  0.  The  Karhunen-Loeve 
decomposition  is  applied  to  the  correlation  function  to  reduce  the  dimensionality  in  the  random  space;  the  generalized 
polynomial  chaos  expansion  is  then  applied  to  the  solution. 

4.2.1  The  Bessel  Correlation  Function 

The  most  commonly  used  correlation  function  for  stochastic  processes  is  the  exponential  function.  In  the  one¬ 
dimensional  case,  it  takes  the  form 

=  (34) 

where  b  is  the  correlation  length.  This  correlation  function  is  the  result  of  first-order  autoregression 

^t  =  a^t-i+et,  a  >  0,  (35) 

where  is  the  random  series  at  t  =  •  •  •  ,  —2,  —1, 0, 1, 2,  •  •  •  and  Ct  is  an  independent  identically  distributed  random 
series.  This  is  a  unilateral  type  of  scheme  where  the  dependence  is  extended  only  in  one  direction,  and  it  is  the 
simplest  realistic  time  series.  For  space  series,  a  bilateral  autoregression  is  more  realistic 

it  =  ait-i  +  bit+i  +  etj  (36) 

where  it  is  intuitively  clear  that  a  and  b  cannot  be  too  large.  It  is  shown  that  the  bilateral  type  of  scheme  is  not 
necessary  in  one  dimension  as  it  can  be  effectively  reduced  to  a  unilateral  one  [35].  Thus  the  exponential  correlation 
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function  can  be  considered  as  the  ‘elementary’  correlation  in  one  dimension.  It  has  been  used  extensively  in  the 
literature  and  its  Karhunen-Loeve  decomposition  can  be  solved  analytically  [11]. 

In  two  dimensions,  the  exponential  correlation  function  can  be  written  as  C(r)  =  where  r  is  the  distance 

between  two  spatial  points.  This  function  has  been  also  used  in  the  literature.  However,  as  Whittle  pointed  out  in 
[35],  it  is  necessary  to  introduce  autoregression  schemes  with  dependence  in  all  directions  for  more  realistic  models 
of  random  series  in  space.  The  simplest  such  model  is 

Cst  =  a(Cs+l,i  +  +  Sst,  (37) 

where  ^st  is  random  field  at  grid  (s,  t)  and  is  independent  identically  distributed  random  field.  This  model 
corresponds  to  a  stochastic  Laplace  equation  in  the  continuous  case: 

^{x,y)  =  e{x,y),  (38) 

where  =  1/a  — 4.  The  ‘elementary’  correlation  function  in  two  dimensions  can  be  solved  from  the  above  equation: 

C(.-)  =  (9 .  (39) 

where  Ki  is  the  modified  Bessel  function  of  the  second  kind  with  order  1,  b  scales  as  the  correlation  length  and  r 
is  the  distance  between  two  points.  On  the  other  hand,  the  exponential  correlation  function  C{r)  =  in  two 

dimensions  corresponds  to  a  rather  artificial  system 

3 

4 

^{x,y)  =  e{x,y).  (40) 

It  is  difficult  to  visualize  a  physical  mechanism  which  would  lead  to  such  a  relation.  For  a  detailed  discussion  on  this 
subject,  see  [35]. 

In  this  paper,  we  employ  (39)  as  the  correlation  function  of  k  and  /.  Since  no  analytical  solution  is  available 
for  the  eigenvalue  problem  (17)  of  the  Karhunen-Loeve  decomposition  for  this  correlation  function,  a  numerical 
eigenvalue  solver  is  employed.  Figure  6  shows  the  distribution  of  the  first  twenty  eigenvalues.  Here  the  parameter  b 
is  set  to  6  =  20.  In  figure  7  and  8  the  eigenfunctions  corresponding  to  the  first  four  eigenvalues  are  plotted. 

4.2.2  Legendre- chaos  and  Uniform  Distribution 

In  this  section  we  assume  k{x,  y;  lo)  and  f{x,  y;  to)  are  random  fields  resulted  from  the  Karhunen-Loeve  decomposition 
(16)  of  the  Bessel  correlation  function  (39),  and  with  the  underlying  random  variables  having  uniform  distributions. 
For  computational  simplicity,  we  further  assume  k  and  /  are  fully  cross-correlated.  Due  to  the  fast  decay  of  eigenvalues 


/ay  /  dV  1 
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Index 

Figure  6:  Eigenvalues  of  KL  decomposition  with  Bessel  correlation  function  (39),  b  =  20. 


Figure  7:  Eigenfunctions  of  the  KL  decomposition  with  the  Bessel  correlation  function  (39),  b  =  20;  Left:  first 
eigenfunction,  Right:  second  eigenfunction.  (Dashed  lines  denote  negative  values.) 


Level  v4 
6  0.0076245 


Figure  8:  Eigenfunctions  of  the  KL  decomposition  with  the  Bessel  correlation  function  (39),  b  =  20;  Left:  third 
eigenfunction,  Right:  fourth  eigenfunction.  (Dashed  lines  denote  negative  values.) 
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as  shown  in  figure  6,  we  choose  the  first  four  eigenmodes  from  the  Karhunen-Loeve  decomposition.  This  results  in 
a  four-dimensional  (in  random  space)  chaos  expansion.  The  corresponding  chaos  in  this  case  is  the  Legendre-chaos 
(table  1). 

The  spectral//ip  element  method  is  used  for  spatial  discretization.  For  detailed  account  of  spectral/ft.p  element 
method,  see  [33].  Specifically,  an  array  of  5  x  5  elements  are  used  in  the  domain  and  sixth-order  polynomials  are 
employed  as  the  (spatial)  expansion  basis  in  each  element  Numerical  tests  show  that  this  is  sufficient  to  resolve 
the  solution  in  space.  The  standard  deviations  of  the  random  inputs  are  1^  =  <^f  =  0.4.  Resolution  checks  in  random 
space  were  conducted,  and  it  was  shown  that  third-order  (p  =  3)  Legendre-chaos  results  in  converged  solution.  For 
4-dimensional  chaos  (n  =  4),  the  total  number  of  expansion  terms  is  35  (see  equation  (23)). 

Since  no  analytical  solution  is  available,  we  employ  Monte  Carlo  simulations  to  validate  the  chaos  solution. 
Here  we  conduct  the  Monte  Carlo  computation  after  the  Karhunen-Loeve  decomposition,  i.e.,  we  generate  the 
random  number  ensemble  on  the  reduced  basis  from  the  Karhunen-Loeve  decomposition.  In  this  way  the  error  from 
generalized  polynomial  chaos  expansion  is  isolated,  while  the  error  introduced  by  the  finite-term  truncation  of  KL 
decomposition,  which  is  well-understood,  is  excluded. 

The  solution  profile  along  the  horizontal  centerline  through  the  domain  is  considered  in  figure  9.  The  mean 
solution  of  Legendre-chaos  and  Monte  Carlo  simulation  with  different  number  of  realizations  are  shown,  together 
with  the  corresponding  deterministic  solution.  A  noticeable  difference  between  the  stochastic  mean  profile  and  the 
deterministic  profile  is  observed.  In  figure  10  the  variance  of  the  stochastic  solution  along  the  horizontal  centerline 
is  shown.  It  is  seen  that  the  Monte  Carlo  solution  converges  to  the  chaos  solution  as  the  number  of  realizations 
increases.  Good  agreement  is  obtained  with  50,  000  realizations. 


Figure  9:  Two-dimensional  model  problem:  uniform  random  distribution  and  Legendre-chaos;  Left:  Mean  solution 
along  the  horizontal  centerline.  Right:  Close-up  view. 

^We  note  here  that  the  Jacobi  polynomials  of  mixed  weights  are  used  in  the  spectral/hp  element  method. 
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Figure  10:  Two-dimensional  model  problem:  uniform  random  distribution  and  Legendre-chaos;  Left:  Variance  along 
the  horizontal  centerline,  Right:  Close-up  view. 


Similar  results  are  obtained  for  other  solution  profiles  in  the  domain,  for  example,  the  vertical  centerline. 

4.2.3  Hermite-chaos  and  Gaussian  Distribution 

We  now  assume  the  random  field  K{x,y;uj)  and  f{x,y;uj)  are  Gaussian  processes  with  =  <Jf  =  0.2.  All  the 
remaining  parameters  are  the  same  as  the  above  example.  The  corresponding  generalized  polynomial  chaos  is  the 
Hermite-chaos . 

The  same  solution  profiles  along  the  horizontal  centerline  of  the  domain  are  shown  in  figure  11  and  12,  for  the 
mean  solution  and  the  variance,  respectively.  In  this  case,  a  fourth-order  Hermite-chaos  (p  =  4)  is  required  to  obtain 
converged  result  in  random  space.  This  corresponds  to  a  70-term  expansion  from  formula  (23)  for  n  =  4,p  =  4.  The 
corresponding  solution  of  the  Monte  Carlo  simulation  converges  relatively  fast  in  this  case,  and  for  20,  000  realizations 
it  converges  to  the  Hermite-chaos  solution. 


Figure  11:  Two-dimensional  model  problem:  Gaussian  random  distribution  and  Hermite-chaos;  Left:  Mean  solution 
along  the  horizontal  centerline.  Right:  Close-up  view. 
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Figure  12:  Two-dimensional  model  problem:  Gaussian  random  distribution  and  Hermite-chaos;  Left:  Variance  along 
the  horizontal  centerline,  Right:  Close-up  view. 

4.2.4  Charlier-chaos  and  Poisson  Distribution 

As  an  example  of  the  discretely  distributed  random  fields,  we  now  assume  the  diffusivity  «:(a;,  y,  lo)  and  source  term 
f{x,y;Lu)  are  processes  resulted  from  Poisson  random  variables  in  the  Karhunen-Loeve  decomposition  (16),  with 
cTk  =  o'f  =  0.2.  The  parameter  A  =  1  as  in  equation  (31). 

The  third-order  (p  =  3)  corresponding  generalized  chaos,  the  Charlier-chaos,  results  in  resolution-independent 
solution  in  random  space.  The  Monte  Carlo  solution  converges  to  the  solution  of  Charlier-chaos;  with  100,000 
realizations  we  obtain  good  agreement.  The  solution  profiles  of  the  mean  and  variance  along  the  horizontal  centerline 
are  shown  in  figure  13  and  14,  respectively. 


Figure  13:  Two-dimensional  model  problem:  Poisson  random  distribution  and  Charlier-chaos;  Left:  Mean  solution 
along  the  horizontal  centerline.  Right:  Close-up  view. 
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Figure  14:  Two-dimensional  model  problem:  Poisson  random  distribution  and  Charlier-chaos;  Left:  Variance  along 
the  horizontal  centerline,  Right:  Close-up  view. 

4.2.5  Krawtchouk- chaos  and  Binomial  Distribution 

Finally,  the  random  field  of  K{x,y;u!)  and  are  assumed  to  have  the  binomial  distributed  random  variables 

with  {N  =  5,  <7  =  0.5)  from  equation  (32)  in  their  Karhunen-Loeve  expansion.  The  standard  deviations  are  = 
af  =  0.2. 

Figure  15  shows  the  mean  solution  along  the  horizontal  centerline  of  the  domain,  while  figure  16  shows  the 
variance  profile.  The  third-order  (p  =  3)  Krawtchouk-chaos  is  sufficient  to  resolve  the  problem  in  random  space.  On 
the  other  hand,  the  solution  of  Monte  Carlo  simulation  converges  to  the  chaos  solution  with  50,  000  realizations. 


Figure  15:  Two-dimensional  model  problem:  binomial  random  distribution  and  Krawtchouk-chaos;  Left:  Mean 
solution  along  the  horizontal  centerline.  Right:  Close-up  view. 
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Figure  16:  Two-dimensional  model  problem:  binomial  random  distribution  and  Krawtchouk-chaos;  Left:  Variance 
along  the  horizontal  centerline,  Right:  Close-up  view. 

4.3  Random  Heat  Conduction  in  a  Grooved  Channel 

In  this  section  we  consider  the  steady  state  heat  conduction  in  a  grooved  channel  subject  to  uncertainties  in  boundary 
conditions  and  diffusivity. 

W  ■[K{x,y;uj)Wu{x,y;uj)]  =  0,  {x,y)  G  D,  (41) 

where  the  computational  domain  D  is  shown  in  figure  17.  The  boundary  of  the  domain  consists  of  four  segments:  the 
top  of  the  channel  Fr,  the  bottom  of  the  channel  Tb,  the  two  sides  of  the  channel  Fg  and  the  boundaries  of  the  cavity 
Fc.  The  diffusivity  K{x,y;uj)  is  a  random  field  with  uniformly  distributed  random  variables  in  its  Karhunen-Loeve 
decomposition,  with  mean  field  K{x,y;uj)  =  1  and  the  same  Bessel  correlation  function  as  in  section  4.2.1.  The 
boundary  conditions  are 

Ou 

w|r.j,  =  0)  ^  ■w|r^=l  +  C,  (42) 

where  ^  is  a  random  variable  with  uniform  distribution.  For  the  spectral//ip  element  solver  in  space,  four  elements 
are  used  in  the  domain,  as  shown  in  figure  17.  Within  each  mesh,  10*^-order  (Jacobi)  polynomials  are  employed. 
In  the  random  space,  the  third-order  Legendre-chaos,  corresponding  to  the  uniformly  distributed  random  inputs,  is 
used.  Resolution  checks  indicate  that  the  above  discretization  is  sufficient  to  resolve  the  problem,  both  in  physical 
and  random  spaces. 

We  consider  two  cases:  the  first  case  is  when  only  the  diffusivity  k  is  random,  while  the  boundary  condition 
along  Fc  is  deterministic,  i.e.  u|rc  =  1-  Same  as  in  section  4.2.1,  the  first  four  eigenmodes  of  the  Karhunen-Loeve 
decomposition  are  employed  to  represent  k.  This  results  in  a  four-dimensional  (n  =  4)  chaos  expansion.  For  third- 
order  chaos  (p  =  3),  a  total  of  35  expansion  terms  are  needed  from  (23).  In  the  second  case,  we  further  assume  the 
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Figure  17:  Schematic  of  the  domain  of  the  grooved  channel 

boundary  condition  along  the  wall  of  cavity  is  random  as  in  (42),  and  is  independent  of  the  random  field  n.  This 
introduces  one  more  dimension  in  the  random  space  and  a  total  of  56  expansion  terms  are  needed  for  third-order 
chaos  expansion;  n  =  5,  p  =  3  from  (23). 

In  figure  18,  the  contours  of  the  standard  deviations  of  the  solution  are  plotted.  The  solution  of  the  first  case 
is  shown  on  the  left,  while  solution  of  the  second  case  on  the  right.  In  both  cases,  the  standard  deviations  of  the 
random  inputs  are  a  =  0.2.  No  noticeable  difference  is  observed  between  the  mean  solutions  of  the  two  cases,  and 
that  of  the  corresponding  deterministic  case.  However,  the  standard  deviations  of  the  solutions  are  very  different  for 
the  two  cases.  From  figure  18,  we  can  see  that  the  effect  of  uncertainty  in  the  diffusivity  is  subdominant  (maximum 
deviation  about  only  0.15%).  By  introducing  the  uncertainty  in  boundary  condition  along  the  walls  of  the  cavity, 
the  output  uncertainty  is  greatly  enhanced  in  the  entire  domain  (maximum  deviation  about  12%),  and  its  structure 
is  changed;  the  maximum  of  the  output  uncertainty  moves  from  the  center  of  the  channel  to  the  lower  wall  of  the 
cavity. 

5  Summary 

We  have  developed  a  stochastic  spectral  method  to  model  uncertainty  in  steady  state  diffusion  problems.  The 
generalized  polynomial  chaos  we  introduced  includes  the  original  polynomial  chaos,  the  Hermite-chaos,  as  a  subset, 
and  is  an  extension  of  the  original  chaos  idea  of  Wiener  (1938)  and  of  the  work  of  Ghanem  &  Spanos  (1991). 
The  important  feature  of  the  new  broader  framework  is  that  it  incorporates  different  types  of  chaos  expansion 
corresponding  to  several  important  distribution  functions,  including  some  discrete  distributions  which  cannot  be 
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Figure  18:  Standard  deviations  of  heat  conduction  in  the  grooved  channel;  Left:  solution  subject  to  random  diffusivity 
only;  Right:  solution  subject  to  random  diffusivity  and  random  boundary  conditions. 


readily  handled  by  the  original  polynomial  chaos  directly. 

We  have  applied  the  generalized  polynomial  chaos  to  the  solution  of  steady  state  random  diffusion  problems, 
as  a  natural  extension  of  our  earlier  work  [17,  19,  20].  In  particular,  we  employed  a  block  Gauss-Seidel  iteration 
technique  to  solve  the  system  of  equations  efficiently.  An  ‘elementary’  correlation  function  in  two  dimensions,  the 
Bessel  correlation  function,  was  studied  and  applied  in  the  computations.  The  Karhunen-Loeve  decomposition  is 
used  to  reduce  the  dimensionality  of  the  random  space. 

We  have  shown  that,  when  the  appropriate  chaos  expansion  is  chosen  according  the  random  input,  the  generalized 
polynomial  chaos  solution  converges  exponential  fast  due  to  the  fact  that  it  is  a  spectral  expansion  in  the  random 
space.  The  exponential  convergence  rate  is  demonstrated  for  one-dimensional  model  problem,  and  is  in  accordance 
with  the  result  of  [17].  For  more  complicated  problems,  the  Monte  Carlo  simulation  is  employed  to  validate  the  chaos 
solution.  We  observe  good  agreement  between  the  well-resolved  chaos  expansion  solution  and  the  converged  Monte 
Carlo  simulation  results.  For  this  particular  problem,  tens  of  thousands  realizations  are  needed  for  the  Monte  Carlo 
computation,  and  the  generalized  polynomial  chaos  expansion  is  at  least  two  to  three  orders  faster.  In  particular, 
we  confine  our  applications  to  problems  within  finite  domain  in  the  physical  space.  Fast  convergence  of  the  chaos 
expansion  is  observed  with  relatively  large  variance  of  the  random  input.  Problems  with  infinite  physical  domain 
require  separate  treatment. 

The  efficiency  of  the  chaos  expansion  is  problem  specific  and  depends  greatly  upon  the  dimensionality  of  the 
random  space.  Although  Karhunen-Loeve  decomposition,  among  other  possible  techniques,  can  be  used  to  reduce 
the  dimensionality,  it  can  be  large  for  systems  with  very  short  correlation  length,  e.g.,  the  white  noise.  To  this  end, 
the  number  of  expansion  terms  may  be  very  large,  thus  reducing  the  efficiency  of  the  chaos  expansion  drasticly.  This 
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problem  also  deserves  further  research. 
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